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Abstract 

We solve the Gross-Pitaevskii equations for a dilute atomic gas in a magnetic 

trap, modeled by an anisotropic harmonic potential. We evaluate the wave 

function and the energy of the Bose Einstein condensate as a function of 

the particle number, both for positive and negative scattering length. The 

results for the transverse and vertical size of the cloud of atoms, as well 

as for the kinetic and potential energy per particle, are compared with the 

predictions of approximated models. We also compare the aspect ratio of 

the velocity distribution with first experimental estimates available for ^^Rb. 

Vortex states are considered and the critical angular velocity for production 

of vortices is calculated. We show that the presence of vortices significantly 

increases the stability of the condensate in the case of attractive interactions. 
PACS numbers: 03.75.Fi, 05.30.Jp, 32.80.Pj 
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I. INTRODUCTION 



Recent experiments support the existence of a Bose Einstein condensate in vapors of 
alkali atoms, such as rubidium litium and sodium 0, confined in magnetic traps and 
cooled down to temperature of the order of 100 nanokelvin. This important discovery opens 
new interesting perspectives in the field of many-body physics (for a comprehensive review 



The vapors of alkali atoms used in the experiments are very dilute, i.e., the average 
distance among the atoms is much larger than the range of the interaction. So the physics 
is expected to be dominated by two-body collisions, well accounted for by the knowledge 
of the s-wave scattering length. This also implies that the Gross-Pitaevskii theory ^ for 
weakly interacting bosons finds in these systems an ideal field of application. This theory 
has been already used to describe bosons confined in isotropic traps The anisotropic 

case has been recently discussed by Baym and Pethick who have obtained approximate 
analytic solutions, thereby providing a first qualitative insight into many interesting features 
of these systems. To make the analysis more quantitative one has to solve numerically the 
Gross-Pitaevskii equations. 

In the present work we provide a complete set of solutions of the Gross-Pitaevskii equa- 
tions applied to alkali atoms in anisotropic traps. We calculate the condensate wave 
function at T = for bosons interacting through positive and negative scattering length. 
We explore the iV-dependence of relevant quantities, such as energy, chemical potential and 
the aspect ratio of the velocity distribution. We also discuss vortex states, studying the 
density profile and the critical angular velocity. 

The results of the numerical minimization of the energy functional are compared with 
the analytic solution of the noninteracting anisotropic harmonic oscillator as well as with 
approximated models available when the interaction is strongly repulsive (strongly repulsive 
limit). For instance we will show that the surface structure of the condensate wave function, 
for which a sound treatment of the kinetic energy is needed, is relevant in determining the 
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aspect ratio. Our results for the aspect ratio in ^ Rb, with of the order of a few thousands, 
are in agreement with the first experimental findings of Ref. An accurate determination 
of the particle distribution in the vapor cloud is also relevant in order to calculate the 
moment of inertia of the system, which is directly connected with the superfluid behavior 

We find instabilities in the solutions for negative scattering length when the number of 
atoms exceeds certain critical values, of the order of 1400 for ''Li, in agreement with previous 
analysis 0. However we find also that the stability of these systems is significantly increased 
when vortex states are considered. This happens because the vortex flow pushes the atoms 
away from the center of the trap decreasing the highest value of the local density. For ^Li we 
find minima of the Gross-Pitaevskii functional corresponding to vortex states with quantum 
of circulation higher than 1 and having of the order of 10000. 

The paper is organized as follows. In Sec. we present the formalism of the Gross- 
Pitaevskii theory for anisotropic traps and we discuss the solutions in two limits (nonin- 
teracting and strongly repulsive limit). We also discuss the generalization of the theory to 
include vortex states. In Sec. Q we briefly introduce the numerical procedure, which is 
based on the steepest descent method for functional minimization. In Sec. we present 
the results for the two cases of positive (^''Rb) and negative (^Li) scattering length. Finally, 
in Sec. we summarize the main results. 



II. GROSS-PITAEVSKII THEORY FOR TRAPPED BOSONS 



In the Gross-Pitaevskii theory the ground state energy for condensed bosons of mass 
m is given by the functional 



E[il)] = Jdr 



where il^{r) is the condensate wave function (order parameter), u!± and are the two angular 
frequencies associated with the external potential of the anisotropic trap and a is the s-wave 
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scattering length. The wave function is assumed to be normahzed to the number of atoms 
in the condensate: 

j dv = N . (2) 

In the T — i> hmit considered in this work, coincides with the total number of atoms 
in the trap. The explicit form of the ground state wave function is obtained by minimizing 
the energy functional. One can also write explicitly the variation of the energy functional 
at first order, finding the Euler-Lagrange equation 



^(r) = /x^(r) , (3) 



where /i is the chemical potential. This equation has the form of a nonlinear stationary 
Schrodinger equation. It has been recently solved for bosons in isotropic harmonic traps 
by Edwards and Burnett 0, by means of iterated Runge-Kutta integrations. An efficient 
numerical calculation of the ground state, which works well also in the case of anisotropic 
traps with and without vortices, consists in minimizing directly the energy functional (|Tp 
with a steepest descent method, as we will show in sect. 

The Gross-Pitaevskii theory is expected to be accurate when the system is dilute. If p is 
the density of bosons, the parameter which measures the applicability of the theory is the 
product a^p, which should be much less then 1. This condition is largely satisfied by the 
samples of alkali atoms used in the recent experiments |l|,0]. For instance, the central density 
of 10000 atoms of ^''Rb in the trap of Ref. [|l| is expected to be of the order of 10^^ -r- 10^^ 
cm~^, which yields a^p ~ 10~^ or less. Even in the experiment of Ref. 0] on sodium, where 
the number of atoms in the condensate is much larger, the quantity a^p remains very small 
(~ 10~^). The theory is also accurate for dilute systems of bosons having negative scattering 
length. In this case, however, one has to take care of the possible instability induced by the 
attractive interaction when is large. When the system collapses, a more accurate theory 
is required in order to include short range effects. 

To simplify the formalism one can choose a slightly different notation, taking advantage 



of the fact that all distances and energies in the calculation scale as the typical length and 
energy of the harmonic external potential. So, we introduce the standard lengths 

a± = (^^] ; az = , (4) 

and we rescale the spatial coordinate, the energy and the wave function as follows: 

r = a^ri (5) 

E = huxEi (6) 



^(r) = yiVM^i(ri). (7) 
The wave function ifji is normalized to 1. Finally we introduce the asymmetry parameter 

A = ujz/^± (8) 

and define the quantity 

SvraA^ 

ui = . 9 



El 

N 



(10) 



With these changes, the Gross-Pitaevskii energy functional takes the form 

and the nonlinear Schrodinger equation becomes 

-V? + xl + yl + Wl + Mi|^i(ri)p] V^i(ri) = 2/iiV^i(ri) , (11) 

where /i = hujj_fii. The dimensionless quantity Ui characterizes the effect of the interaction 
in the equation for the condensate. It is worth noting that even if the system is very dilute 
{a^p <^ 1) the interactions can nevertheless play a crucial role in determining the solution 
of the Gross-Pitaevskii equations. In fact the two conditions a^p ^ 1 and mi ~ 1 can be 
well satisfied for realistic values of the parameters of the problem, i.e., the scattering length, 
the oscillator frequencies and the number of atoms in the trap. We now discuss briefly the 
two limiting cases of noninteracting particles and of strongly repulsive interactions, where 
the solution of Eq.(|TT]) is available in analytic way. We then show how the formalism can 
be extended to describe vortex states. 



A. The noninteracting model 



When the scattering length a vanishes, the problem reduces to the solution of the sta- 
tionary Schrodinger equation for an anisotropic harmonic oscillator. The ground state wave 
function is 



exp 



(12) 



The chemical potential coincides with the energy per particle, which turns out to be (H-A/2). 
The gaussian has different transverse and vertical widths. In particular one has (xf) = (yf) = 
1/2 and (zf) = 1/(2A). The mean values of the momentum operators and pi can be also 
easily calculated. In particular it is interesting to calculate the quantity 

^{Pl)/{PI) = = (13) 



In the following we consider the quantity y {pD / {pi) as a measure of the aspect ratio which 
characterizes the anisotropy of the velocity distribution. This is a relevant quantity in the 
interpretation of the experimental results. Values of the aspect ratio different from 1 reflect 
a peculiar and unique feature of Bose Einstein condensation. 



B. Strongly repulsive limit 

The opposite limit is obtained when the interaction is so strong, or the number of particles 
so large, that the kinetic energy can be neglected in the energy functional. It corresponds to 
very large values of the parameter ui (see Eq. (^). This limit has been already discussed in 
Refs. HJ^. The solution is easily obtained by dropping the kinetic energy term in Eq. (|ll]). 
It has the form 

i,l{r{) = -{2^,-xl-yl-\'zl) (14) 

Ui 

if the right hand side is positive, and ipi = Q elsewhere. The chemical potential is easily 
calculated by imposing the normalization condition / ipfdri = 1. One finds 
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/ 1 ^ \ 2/5 

2/ii = (^Awij , (15) 

and /i = hu±fii. Using the definition of Ui given in Eq. (P) and the relation = dE/dN, 
one also gets the relationship E/N = (5/7) /i. The cloud of atoms extends over a radius 
Ri = y/2fii and a vertical size Zi = XRi. One also finds the results 

i^l) = ^ ; (4) = ^ (16) 

for the square radius (xf) and (zf). Due to the different scaling properties of the wave 
function with respect to the variable z (compare Eqs. ([T2|) and (0)), the aspect ratio 
\/ {pD I kVx) this case is equal to A and not to \f\ as in the noninteracting case. The 
central density of the cloud is ^2yUi/Mi. The wave function (|T3|) is expected to approximate 
well the exact solution of the nonlinear Schrodinger equation ([n|) for large A^, apart from 
the structure of the surface region where the exact wave function has to vanish smoothly. 
This fact has been already tested for isotropic traps in Ref. ||^. However some relevant 
observables can be significantly affected by this surface structure even at large A^, as we will 
see later. 



C. Vortex states 

The energy functional ([T0|) is easily generahzed to include vortex states, i.e., states 
where all the atoms rotate around the z-axis with quantized circulation. Indeed one of the 
primary motivations of the Gross-Pitaevskii theory was the study of vortex states in weakly 
interacting bosons . One has to consider a complex condensate wave function of the form 

^(r) = V^(r)exp[i5(r)] (17) 

where ^/'(r) = y^p(r) is the modulus, while the phase S acts as a velocity potential: v = 
{h/m)'WS. By choosing S = K,(f), where (f) is the angle around the z-axis and k is an integer, 
one has vortex states with tangential velocity 
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with = + y'^. The number k is the quantum of circulation and the angular momentum 
along z is Nnh. Now, one has to put the complex wave function \E' in place of in the energy 
functional (|I]). Using the adimensional quantities defined in Eqs. (|^, |^) and rf^ = xl + y\, 



the resulting Gross-Pitaevskii functional becomes 
dri 



El 

N 



(19) 



which differs from functional (|T0|) because of the centrifugal term. The corresponding non- 
linear Schrodinger equation is 



VI + + rt± + X'zt + MilV-iln)! ^i(ri) = 2/ii^i(ri) . 



(20) 



Due to the presence of the centrifugal term, the solution of this equation for k, ^ has to 
vanish on the 2;-axis. 

For noninteracting particles one falls again in the case of the stationary Schrodinger 
equation for the anisotropic harmonic potential. For instance the k = 1 solution has the 
form 



?/'i(ri) cx ri_Lexp 



(21) 



The energy per particle of the k, ^ states of the anisotropic harmonic oscillator is simply 
Khu± plus the ground state energy. In our dimensionless notation: /xi = 1 + (A/2) + k. 

In the interacting case the kinetic energy can not be neglegted even for large N, since it 
determines the structure of the vortex core. In particular, the balance between the kinetic 
energy and the interaction energy fixes a typical distance over which the condensate wave 
function can heal. For a dilute Bose gas the healing length is given by H 



^ = (Svrpa) 



-1/2 



(22) 



where p is the density of the system. In the case of a vortex it corresponds to the distance 
over which the wave function increases from zero, on the vortex axis, to the bulk density. 
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For the trapped atoms in the N ^ oo hmit we have seen that the central density of the 



cloud is about y2/ii/Mi, where ui and /ii are given in Eqs. @ and (0), so that the healing 
length (in units a±) is ~ (2/ii)~^/^. Since the radius Ri of the cloud, in the same units, 
is of the order of one has 

£i 1 1 



Ri 2/ii Rf 

or, equivalently. 



(23) 



R \RJ 



(24) 



Thus the healing length is small compared with the size of the cloud if R is much bigger 
than a±. 

Vortex states play an important role in characterizing the superfluid properties of Bose 
systems, as is well known in the case of superfluid Helium [|TT|]. The critical angular velocity 
required to produce vortex states is easily calculated once the energies of the states with and 
without vortices is known. One has to compare the energy of a vortex state in frame rotating 
with angular frequency Q, that is {E — QL^), with the energy of the ground state with no 
vortices. Since the angular momentum per particle is nh, the critical angular velocity, in uj± 
units, is 

= [{Ei/N)^ - {Ei/N)o] . (25) 

In the noninteracting case the difference of energy per particle is simply k, so that fic = 1; 
the critical angular velocity is just the angular frequency of the trap in the x, y plane. 



III. NUMERICAL PROCEDURE 

The main purpose of this work is the numerical minimization of the Gross-Pitaevskii 
functional (p!9| ) in order to calculate the properties of the ground state (k = 0) and of 
vortex states (k 7^ 0) for given values of the parameters A^, A, a±_ and a. A method of direct 
minimization is provided by the steepest descent approach. In brief, it consists of projecting 
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onto the minimum of the functional an initial trial state by propagating it in imaginary time. 
A time dependent wave function ipi{ri,t), where t is a fictitious time variable, is evaluated 
at different time steps, starting from an arbitrary trial function and converging to the exact 
solution ipiiri, oo) = ipi{ri). The time evolution can be formulated in terms of the equation 

dt 5?/^i(ri,t) 

where S indicates the constrained functional derivative that preserves the normalization. 
This equation defines a trajectory in the wave function space (and the fictitious time is 
just a label for different configurations) in which at each step one moves a little bit down 
the gradient —EE/Sip. The constrained functional derivative is obtained by adding the 
normalization condition to the functional derivative 

''^'/'^ -HMr,.t), (27) 



where H depends nonlinearly on ipi. The end product is the self-consistent minimization 
of the energy, which corresponds to (dipi/dt) = or, including the normalization, to the 
equation Hipi = 2^iipi, which coincides with the nonlinear Schrodinger equation (^^. In 
practice one chooses an arbitrary time step At and iterates the equation 

^i(ri, t + At) ~ ^i(ri, t) - At Ht/j^iru t) , (28) 

by normalizing ipi to 1 at each iteration. The time step At controls the rate of convergence. 



Several methods have been proposed in the recent literature (see for instance Ref. |T2 
and reference therein) to improve the steepest descent method of functional minimization; 
however, the Gross-Pitaevskii functional is much simpler than the typical functionals used 
in strongly correlated systems and the steepest descent method described above is enough 
efficient for our purposes. 

In practice one has to discretize the {ri±, z)-space with a two-dimensional grid of points, 
so that the wave function becomes a matrix. At each time step the matrix elements are 
changed as in Eq. (^), where the derivatives entering the Hamiltonian are evaluated by 
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means of finite difference formulae. The algorithm can be tested by comparing the results of 
the noninteracting case with the analytical solution of the anisotropic harmonic potential. 
In the interacting case, with large A^, it is also possible to compare the numerical results with 
the analytic solution ([T^). Another test of accuracy is given by the virial theorem, which 
fixes rigorous relationships among the different contributions to the kinetic and potential 
energy of the system at any value of A^. 

The system is sufficiently well described using a grid of 50 x 50 points in the range 
< ri± < 5, and the same for z. The number of iterations in imaginary time depends 
on the degree of convergence required and the goodness of the initial trial wave function. 
The latter can be one of the two analytical limits already discussed, but the final results 
do not depend on the trial wave function. Typically we use 2000 -t- 10000 iterations. Since 
the internal energy is a local functional, each iteration is very fast, so that the functional 
minimization takes no more than 2-^3 minutes of CPU on a DEC- Alpha processor. 



IV. RESULTS 



A. Positive scattering length: Rb 

As an example of atoms with repulsive interaction we choose ^^Rb, as in the experiment 
of Ref. [|T|. The s-wave triplet-spin scattering length is in the range 85ao < a < 140ao, where 
is the Bohr radius DTSf. In our analysis we use a = IOOoq. The asymmetry parameter of 
the experimental trap is A = uJz/uj± = V8. The axial frequency lu^ /2Tr is taken to be 220 Hz 
T^ . The corresponding characteristic length is a± = 1.222 x 10^^ cm and the ratio between 



the scattering and the oscillator lengths is a/a± = 6.47 x 10^^. 

We minimize the Gross-Pitaevskii functional in a wide range of particle number A^. 
Results for the chemical potential and the energy per particle are shown in Table |. Both 
quantities are expressed in units of hu^, or of the equivalent temperature huj±/kB = 3.73 nK. 
The partial contributions to the energy per particle coming from the kinetic energy (kin), 
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the harmonic oscillator potential (ho) and the internal potential energy (pot) are also given. 
The = 1 case coincides with the noninteracting anisotropic harmonic oscillator: in this 
case the internal potential energy vanishes, the kinetic energy and the harmonic oscillator 
potential energy are equal, the chemical potential and the total energy per particle are both 
equal to the analytic value (1 + A/2)=2.41. When increases the repulsion among atoms 
tends to lower the central density, expanding the cloud of atoms towards regions where the 
trapping potential is higher. This produces an increase of both the internal and the harmonic 
oscillator potential energy per particle. Conversely, the kinetic energy per particle decreases 
because the density distribution is flattened. In the strongly repulsive limit, N —>■ oo, one 
should find that the internal potential energy is much greater than the kinetic energy, which 
is the case discussed in Sec. |III3| . Indeed the convergence towards this limit turns out to 
be rather slow. An approximate estimate of the kinetic energy per particle can be obtained 
assuming the wave function to be a gaussian, having a width of the order of the radius R 
of the cloud 0. In this model the kinetic energy is of the order of h'^ / {2mR'^) , which is 
much smaller than the internal potential energy even for relatively small A^ (a factor 10"^ 
is reported in Ref. for A^ = 2000). The discrepancy between the gaussian approximation 
and the exact solution is well understood by looking at the effect of the surface structure 
of the cloud. In Fig. |l] we plot the profiles of the wave function along the x- and 2;-axis for 
several values of A^. The noninteracting case is shown as a dashed line. Increasing A^ the 
central density is significantly lowered. The density in the cloud becomes almost flat and 
it is well approximated by the analytic solution (|1^) valid in the strongly repulsive limit. 
At the surface the wave function vanishes gradually, the typical decay length being almost 
independent of A^. The contribution of the surface to the kinetic energy remains sizable even 
for large A^, so that the kinetic energy is larger than the gaussian estimate h^/{2mR'^). A 
typical profile of the condensate wave function tpi is plotted along the x-axis for A^ = 5000 
in Fig. 1^. The exact minimization of the Gross-Pitaevskii functional (solid line) is compared 
with the noninteracting case (dashed line) and the strongly repulsive limit (dot-dashed). 
Simple relationships among the different contributions to the total energy are obtained 
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by means of the virial theorem . When apphed to the anisotropic trap it gives the rigorous 
relation 

(pi) m , 1 , , 

and analogously for y and z. Summing over the three equations for x, and z one finds 

2-Efcj„ — 2Eho + <^Epot = . (30) 

One can easily see that the numerical results in Table | agree very well with this relation. 

The average size of the cloud in both directions can be easily evaluated once the ground 
state wave function is known. In the last two columns of Table | we report the quantities 



{xf) and Y (zf). When increases the quantity (xf) deviates rapidly from the noninteract- 
ing value 1/2, reflecting the spreading of the atom distribution in the direction of the softer 
trapping potential. The increase of (zf) is slower, but never negligible. On can compare the 
results of the numerical solution of the Gross-Pitaevskii equations with the ones obtained 
in the strongly repulsive limit (see Eq. (0)). In Table we give the approximated chemi- 
cal potential ([T5| ) and the average sizes (0), using the same input parameters (frequencies 
of the trap and scattering length). Comparing these values with the ones in Table |I|, one 
clearly sees that the strongly repulsive limit provides good estimates for the quantities {xD 
and (zf). This means that the behavior of the surface structure, which is very different in 
the exact and approximated wave functions, does not affect significantly the average sizes of 
the cloud. Actually, the estimate of (xf) is better than the one for (zf), since the exact wave 
function approaches more rapidly the one of the strongly repulsive limit in the direction of 
the softer trapping potential. The approximated values of the chemical potential are close 
to the exact ones for N very large. The quality of the strongly repulsive approximation is 
improved in systems with greater values of the parameter ui, as in the case of the sodium 
vapor used in the experiment of Ref. 0], where ~ 10^ and Ui is of the order of 10^. 

Another interesting quantity which can be easily calculated from the ground state wave 



function is the aspect ratio of the velocity distribution, that is the ratio y{pl)/{pl)- This 
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quantity is equal to \/~X in the noninteracting case and should approach A in the strongly 
repulsive limit. The numerical results, as a function of A^, are shown in Fig. |^. The two 
limiting cases are shown as dashed lines. One clearly sees that the convergence to the value 
2.828 = A is very slow; the aspect ratio remains well below the asymptotic value even for 

= 20000. The aspect ratio measured in Ref. is estimated to be about 50 % larger 
than the noninteracting value, while the number of particles is of the order of 5000. The 
agreement with our results is good, even if one has to consider that the experimental estimate 
implicitly assumes a ballistic expansion of the atoms after switching off the external trap. 
The effects of the interaction on the expansion of the gas should be explicitly taken into 
account in order to draw more definitive conclusions. 

Let us now consider the vortex states. In Fig. ^ we show the wave function of a cloud 
of 5000 atoms; the n = 1 wave function (b), which corresponds to atoms flowing around 
the 2;-axis with angular momentum Nh, is compared with the k = ground state (a). The 
atoms are pushed away from the axis forming a toroidal cloud. From the energy of the 
vortex states we calculate the critical angular velocity, through Eq. (pS)). The results for 
n = 1 are shown in Fig. ^. The critical angular velocity decreases rapidly with A^. For 

> 5000 it is less than 40% of the noninteracting value, given by the transverse angular 
frequency uj^ of the trap. A rough estimate of the critical frequency in the large A^ limit 
is given by ^ VLc/uj±_ ~ (a^/_R)^ ln(_R/^), where R is the radius of the cloud. The healing 
length is the distance over which the wave function grows from zero to the bulk value. In 
the limit of large systems it can be approximated by Eq. (|22D with p equal to the density in 
the central part of the toroidal distribution. Both the estimate of ^ and Vic obtained in this 
way are in qualitative agreement with the behavior of the numerical solutions. One can also 
find solutions for k > 1. The critical velocity turns out to increase with k. For instance, 
the critical frequency i7c/2vr for the creation of vortices in a system of 10000 atoms is 26, 35 
and 41 Hz for k = 1, 2 and 3, respectively. 

Finally it is worth recalling that the dimensionless parameter characterizing the effects of 
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the interactions in the Gross-Pitaevskii equations is given by ui = SnaN/ax (see Eqs. (|9|,[T0D ) . 
This imphes that all the results obtained in the present work can be applied, with a proper 
rescaling of the variable A^, to different choices for a and/or a±. For instance, changing the 
axial frequency of the trap from 220 Hz to 120 Hz, so that a± increases by a factor y/ll/6, 
is equivalent to keeping a and a± unchanged and reducing the number of atoms by the same 
factor Jll/Q. 



B. Negative scattering length: "^Li 

As an example of atoms with attractive interaction we choose ^Li, as in the experiment of 
Ref. [@]. The s-wave triplet-spin scattering length is — 27ao ||15|. The axial frequency reported 



in Ref. 0] is uJz/2tt = 117 Hz and the corresponding characteristic length is a± = 2.972 x 10~^ 
cm, thereby yielding a ratio |a|/a_L = 0.48 x 10^'^. The transverse frequency is u}z/2tt = 163 
Hz, so that the asymmetry parameter is A = uJz/ujx_ = 0.72. 

The first important point to stress is that Gross-Pitaevskii functional has no global min- 
imum for negative scattering length. This reflects the tendency of the system to collapse. 
For spatially inhomogeneous systems, however, the zero-point energy can exceed the attrac- 
tive potential, producing local minima of the functional when the density of atoms is not 
too high. The nonlinear stationary Schrodinger equation provides the solutions ip for which 
SE/Stp = 0, but does not say anything about the stability of these solutions. A proper 
treatment of the stability requires a time dependent theory 0J^. The minimization of 
the Gross-Pitaevskii functional with the steepest descent method explores the configuration 
space with axial symmetry near the local minimum. 

In Fig. P we show the results for the wave function along the x- and z-axis for several 
values of A^. As in Fig. ^ we plot the noninteracting case with a dashed line. Here the vapor 
extends more along z than along x, just because the external potential of the trap is softer 
in the axial direction (A < 1). Apart from this purely geometrical fact, the most striking 
difference with respect to the repulsive case is that here the central density of the cloud 
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increases rapidly with A^. This is the effect of adding more and more attractive potential 
energy. When the central density reaches a certain critical limit the system collapses. In 
term of functional minimization this implies that the convergence towards the local minimum 
becomes slower and slower, until a critical above which the energy falls down and does 
not converge anymore. In ^Li, with the input parameters given above, the critical number 
turns out to be about 1400. 

Looking at Fig. ^ one also notices that the wave function changes its form in the same way 
along X and z. Both the average sizes \J (x^) and ^ {z^) decreases slowly when N increases. 
For instance, for = 1000 one has \J {xf) = 0.62 and \f{zf) = 0.69, both values being about 
15% smaller than the ones in the noninteracting case. The ratio ^ {x"^) / {z"^) is practically 
independent of A^. For instance for A^ = 1 it is equal to = 0.85, while for A^ = 1000 
it is 0.90, with an increase of only 5%. The aspect ratio of the velocity distribution, which 
is equal to \/~\ in the noninteracting case, behaves in the same way. Even the energy per 
particle depends smoothly on A^. In units Tioj^, it is equal to 1.36 and 1.15 for A^ = 1 and 
A^ = 1000, respectively. 

Coming back to the question of the stability, we notice that, when the local minimum 
associated with wave functions of the form shown in Fig. |^ disappears, nothing prevents 
a priori the existence of other local minima associated with different configurations. Such 
configurations should have local density lower than the critical one. A natural way to obtain 
a favorable situation is to move the atoms away from the z-axis, conserving the total number 
of particles. This happens in the presence of a vortex. In Fig. ^ we show the wave function 
for 1000 ''Li atoms with no vortices (a) and with an axial vortex of unit circulation (b). 
We use the same units in both cases, so one can see that the maximum value of the wave 
function inside the toroidal distribution of the vortex is approximately a factor two lower 
than the central value in the state with no vorticity (the density is four times smaller). The 
critical angular frequency for the formation of the vortex state in Fig. |^ is 1.12 times the 
transverse angular frequency of the trap. In systems with attractive interaction the critical 
angular velocity is larger than for noninteracting particles, while the opposite is true for 
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repulsive interaction. This is because it costs internal potential energy to lower the average 
density, as the vortex does, for attractive interactions. However, once a vortex is created, 
the corresponding state is more stable than in the absence of vorticity: one can put more 
atoms inside the rotating cloud before reaching the critical density for the final collapse. 
Indeed we find local minima of the Gross-Pitaevskii functional for much larger than 1400 
if K > 0. We show three examples in Fig. |[ One notices that the maximum density of 
the K = 1 state slightly increases from the case = 1000 (Fig. 0b) to the case N = 3500 
(top of Fig. H), however remaining well below the value of the central density in the state 
without vorticity (Fig. |^). The three vortices in Fig. | have almost the same peak density, 
but very different number of particles. They correspond to well defined local minima of the 
functional. If the number of particles is increased, one finds again critical values of for 
which the minima disappear. For k = 1 we find a critical value of ~ 4000; for k, = 2 
and 3 we find critical values of 6500 and 8300, respectively. It is worth mentioning that the 
number of particles in the condensate reported in the experimental work of Ref. [Q] is an 
order of magnitude higher than the critical value for the stability of the Gross-Pitaevskii 
solution without vorticity (A^ ~ 1400). This discrepancy between the experimental finding 
of Ref. and the predictions of the Gross-Pitaevskii theory could be significantly reduced 
if one assumes the existence of a vortex in the atomic cloud. 

V. CONCLUSIONS 

In this paper we have solved the Gross-Pitaevskii equations for a dilute gas of alkali 
atoms in anisotropic magnetic traps by numerical minimization of the total energy. The 
theory provides the condensate wave function at T = for states with and without vorticity. 
The comparison with approximate models, such as the noninteracting gas and the strongly 
repulsive limit, has been carefully explored. We have explicitly discussed the results for ^''Rb 
(positive scattering length) and ''Li (negative scattering length) since these elements have 
been recently used in the first successful measurements of Bose Einstein condensation [|T],0 . 
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We summarize here the main results of the present analysis. 

• We have explored in a systematic way the density distribution and the energy system- 
atics of the atomic clouds. The exact condensate wave function is flat in the interior 
and vanishes smoothly at the surface. The contribution of the surface to the kinetic 
energy per particle remains sizable even for relatively large A^, differently from the 
predictions of approximated models recently proposed. This affects significantly the 
behavior of the aspect ratio of the velocity distribution. In the case of positive scat- 
tering length the aspect ratio is larger than the value a/A given by the noninteracting 
model but smaller than the value A given by the strongly repulsive limit. The values 
calculated for ®^Rb are in agreement with the experimental findings of Ref. 

• We have studied the properties of vortex states. For systems with repulsive interaction 
the critical angular velocity for the formation of vortices decreases rapidly with 
with respect to the value of the noninteracting gas. Conversely, it increases with N 
in systems with attractive interaction. The most striking feature of vortex states is 
the tendency to lower the peak density in the cloud of atoms. This tendency has a 
dramatic effect for systems with attractive interaction, where high values of the peak 
density can produce the disappearance of the local minimum of the functional, i.e, 
the collapse of the system. In ''Li this happens for ~ 1400. It turns out that 
the presence of a vortex increases the stability of the system, in the sense that local 
minima with larger N can be found. We have shown the results up to = 8000 with 
circulation number k = 3. Higher values of A^ can be obtained by increasing k. 

Further work is planned in order to study in more details the velocity distribution of 
the atomic vapor. Time dependent calculations are also feasible within the same theoretical 
scheme. 
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FIGURES 

FIG. 1. Ground state wave function for ®^Rb along the x-axis (upper part) and along the z-axis 
(lower part). Distances are in units a± (see Eq. (^)). The dashed line is the noninteracting case; 
the sohd lines corresponds to = 100, 200, 500, 1000, 2000, 5000 and 10000, in descending order of 
central density. 

FIG. 2. Ground state wave function for 5000 atoms of ^'^Rb. Dashed line: noninteracting case 
(see Sec. [I A). Dot-dashed line: strongly repulsive limit (see Sec. II B). Solid line: exact solution 



of Gross-Pitaevskii functional. 

FIG. 3. Ratio of the axial to transverse average velocity as a function of in ^'^Rb. The lower 
and upper dashed lines corresponds to \/A and A, respectively. 

FIG. 4. Wave function, in arbitrary units, of 5000 ^'^Rb atoms, a) Ground state, b) Vortex 
state with k = 1. 

FIG. 5. Critical angular velocity, in units uj±, for the formation of k = 1 vortices in ®^Rb vapor 
as a function of A^ 

FIG. 6. Ground state wave function for ^Li along the x-axis (upper part) and along the z-axis 
(lower part). Distances are in units a± (see Eq. (§)). The dashed line is the noninteracting case; 
the solid lines corresponds to A^ = 200, 500 and 1000, in ascending order of central density. 

FIG. 7. Wave function, in arbitrary units, of 1000 ''Li atoms, a) Ground state, b) Vortex state 
with K = 1. 

FIG. 8. Vortex state wave functions, in arbitrary units, for different values of A^ and k in 
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TABLES 



TABLE L Results for the ground state of ^"^Rb atoms in a trap with A = Vs. Chemical 
potential and energy in units huj±, with 27ra;_L = 220 Hz. Lengths in units aj_. 



IN 




) 


( rp /]\T\ 




jpot 


1 lrA\ 

V 




1 


Z.4i4 


O A^ A 


1 007 


1 on'7 


n nnn 
U.UUU 


U. / U ( 


n /1 9n 
U.4zU 


iUU 


Z.oo 




i.Uo 


i.oy 


n 91 


n vn 

u. ( y 


n /I /I 

U.44 


onn 

ZUU 


O.Zi 


Z.oD 


u.yo 


1 

i .OZ 


n 

U.oD 


U.oO 


n /I 

U.40 


DUU 


Q/l 


O.oU 


U.oD 


1 SI 
i.oi 


U.Do 


n QR 
U.yo 


n /1 7 

U.4 ( 


1 nnn 

iUUU 


A 77 


0.o4 


U. / 


9 11^ 
Z. 10 


u.yo 


1 ns 

i.Uo 


n i^n 

U.OU 


2000 


5.93 


4.61 


0.66 


2.64 


1.32 


1.23 


0.53 


5000 


8.14 


6.12 


0.54 


3.57 


2.02 


1.47 


0.59 


10000 


10.5 


7.76 


0.45 


4.57 


2.74 


1.69 


0.65 


15000 


12.2 


8.98 


0.41 


5.31 


3.26 


1.84 


0.70 


20000 


13.7 


9.98 


0.38 


5.91 


3.68 


1.94 


0.73 



TABLE IL Chemical potential (in units hLO±^) and average transverse and vertical size (in units 
aj_) in the strongly repulsive approximation (Sec. II B| ), for ^''Rb in the same trap of Table |. 
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